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Abstract 

3D numerical simulations of ferromagnetic materials can be compared with 
experimental results via microwave susceptibility. In this paper, an optimised 
computation of this microwave susceptibility for large meshes is proposed. The 
microwave susceptibility is obtained by linearisation of the Landau and Lifchitz 
equations near equilibrium states and the linear systems to be solved are very ill- 
conditionned. Solutions are computed using the Conjugate Gradient method for 
the Normal equation (CGN Method). An efficient preconditioner is developed con- 
sisting of a projection and an approximation of an "exact" preconditioner in the set 
of circulant matrices. Control of the condition number due to the preconditioning 
and evolution of the singular value decomposition are shown in the results. 

1 Introduction 

Ferromagnetic simulation via the micromagnetic model is a real-life computational 
challenge. Ferromagnetic materials are used in numerous applications such as radar 
protection, magnetic recording or micro electronics. In these applications, the mag- 
netic objects studied are micro or nano-objects which are difficult and expensive to 
craft. Thus, one of the optimisation solutions, for the shape and composition of such 
particles, is numeric simulation. The first step in this type of simulation is to com- 
pute the dynamic of the magnetisation and the equilibrium states. However, a direct 
comparison of the results with experiments is impossible for 3D particles. The main 
comparison tool is microwave susceptibility as the resonnance numerical curves can 
be compared with the physical experiments. At that point several difficulties are en- 
countered. The main one is managing a large number of degrees of freedom. This is 
required to compute interesting configurations with sufficient accuracy. 

In this article, we use the micromagnetism model in order to model the magneti- 
sation behaviour in ferromagnetic materials. This model is a mesoscopic model, ie. a 
model valid for a scale between the one used for microscopic Maxwell equations and 
the scale of classic macroscopic Maxwell equations. In this model, magnetisation does 
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not linearly depend on magnetic excitation but is controlled by a non- linear system: the 
Landau-Lischitz equation 0. This model was introduced by Brown (012]. 

There are two ways to obtain the equilibrium states. The first by energy minimi- 
sation |5J), the second by relaxation of the dynamic system (|6 71)- The main 
advantage of the dynamical approach is to compute an equilibrium state linked to given 
initial data by a life-like dynamic process; then, we can apply dynamical treatments, 
via the external field, in order to find specific equilibrium states. 

Computation of the microwave susceptibility can be performed by two main meth- 
ods: the Harmonic Direct Computation and the Fourier Transform Method. The first 
method is based upon the use of a linearised version of the evolution equation pertu- 
bated by a time harmonic external field. The second is based upon the injection of an 
harmonic perturbation. The Fourier method implies the resolution of a time dependant 
problem that is quite ill-conditionned for low frequencies (the time step ensuring that 
the convergence vanishes swiftly when the frequency decreases) but the linearisation 
methods permit the range of frequencies used in the applications to be attained. 

2 The microwave susceptibility problem 
2.1 The linearisation 

In this problem, we are interested in computing the microwave response of a ferro- 
magnetic system to an external harmonic exitation. We consider that the ferromagnetic 
material is homogeneous and contained in a C 1 -class piecewize domain of R 3 denoted 
SI. Then, we study the evolution of the magnetisation field in the neighbourhood of 
an equilibrium state of the dynamic equation. This equation, in the micromagnetism 
model 1 1 1, is given by the Landau-Lifchitz system: find m in i? 1 ([0, T] x R 3 , SI; R 3 ) 
= {m G L 2 (R 3 ;R 3 )|Vi G [0,T], m\ a G ff^tyR 3 ) and m = in R\Sl} such that 

/ HST ~ /( m > h ^t) = ~ m A (H(m) + £) - a m A (m A (H(m) +£)),€ (0, T] x SI, 
\ m(x,0) = m {x), Vi G SI. 

where H is a linear operator, from i^QO.T] x R 3 ,ft;R 3 ) into ff-^R 3 ; R 3 ), £ the 
external magnetic field (independent of the magnetisation and element of L°°([0, T] x 
R 3 ; R 3 ), a the damping factor (a strictly positive real) and mo is a given element of 
S 2 (Sl)={m G 5" 1 (R 3 ,12;M 3 ) | |m|n| = 1, a.e. in SI}. In this model, we can see that 
the local module of the magnetisation is naturally preserved. In this article, we define 
H as follows: Vm £ F x ([0, T] x R 3 , SI; R 3 ) 

H(rn) = AAm + Hd(m) + K(m — (m.u)u) 

where A and K positive real constants and u is an element of i? 1 ([0, T] x R 3 , SI; S 2 ) 
(S 2 designates the unit sphere). The operator H c i is defined in the sense of distributions 
on R 3 by 

f mt(ffj(m)) = 0, 

[ div(Hd(m)) = — div(m). 

Now, let us define the equilibrium states of the system Q 
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Definition 1 For a given lin L°° (R 3 ; R 3 ) ( independent of time ), a magnetisation state 
mi, in iJ^R 3 , f2; R 3 ) is an equilibrium state if, and only if, 

f(me,£) — 0, a.e. in Q. 

Then, for a given equilibrium state rri£, associated to an external state I, we define the 
microwave susceptibility 

Definition 2 For a given equilibrium state me, associated to an external field I, we 
denote a susceptibility tensor of the order 3 complex matrices defined by 

(xWkk = - 2^(A fc ,mi)o,n, V(l,fc) G {1,2,3} 2 , 

with Afc = C,ke. lut and C,k is a contant vector q/R 3 . Furthermore, we suppose that 
(Cfc)fce{i,2,3} is an orthogonal basis q/R 3 . Then, for all k in {1,2, 3}, m k is a solution 
o/0/or the external field Xk + I and the intial data m — m i- 

Formally, if the excitation is sufficiently small, then the magnetisation responses 
will be also small and we can define this response for every k in {1, 2, 3} by 

m k -m e = pb k e iwt , 

with pk G H 1 (R 3 , fi; C 3 ). In the following we suppose that C,k and p,u are of the same 
order. 

Then, if we re-write the system Q verified by m k , the linearised equation gives 

(iuj - D lt i oh- D 2 j)(p k ) = Di 4 {C,k) (2) 
where, for all w in i°°(R 3 ; R 3 ), we set 

Di,e{w) — —mi A w — a me A (rni A w), 
D 2 ,e(w) = (H(m e ) + £) A w + a m e A (w A (H(m e ) + £)) 



2.2 The discretisation of the linearised equation 

In order to discretise the equation, we consider a monolith K (O) such that ft C K (f2). 
Ideally, this monolith is the smaller containing f2. Then, K(fl) is discretised using a 
regular cubic mesh of cells (Oj)j £ jv h where h is the length of a cell and is the set of 
the indices. We set flh — [J i£N . t h f2j where Ni n t,h C Nu is the set of indices such 
that, for every i in Ni ntt h, f2j fl f2 ^ 0. 

Then, we choose as a discrete space for all euclidian space F: 

W h (F) = {uE L 2 (R 3 ; F)\u = in R 3 \K(n) and Vi G iVfo, U| 0l is a constant}, 

for each u in Wh, we set: Vi G Nh, Ui = u^.. We choose the L 2 scalar product on 
R 3 as the scalar product on Wh, we denote it (u, u)o,n f° r all u,v in L 2 (R 3 ; F). Then, 
setting 

Ph 

L 2 (R 3 ;F) — » (F) 

u .— > i\(u) = ^ (-1 / wdx) 
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where lj is defined for x in R 3 by \i(x) — 1 if x belongs to Vli, li(x) = otherwise. 
P£ designates the canonical injection of Wh{F) onto L 2 (M. 3 ; F). 

These definitions lead to the following formulas for the discrete magnetic contribu- 
tions: 

H a , h =P h oH a oP*, 

and 

H dth = P h oH d oP*, 

the analysis of H a ^h is straightforward. On the other hand, the analysis of H d ,h is 
not direct, in particular, it has been demonstrated that this discretisation preserves the 
main properties of the demagnetisation operator Hd {Hd is a projection operator), and 
a lower estimate of its lower eigenvalue is given. Furthermore, the computation of this 
operator is very expensive: the discrete matrix is a full matrix. Then, to optimise its 
computation, we choose to use a regular cubic mesh which ensure a specific structure 
for the discrete operator. This block-Toeplitz structure enables us to reduce the storage 
of the matrix from #(Nh) 2 to 0(#(JV/,)) and the computation cost from #(Nh) 2 to 
0(#(JV/j) log(#(iV/i))). For complete analysis of the discretisation of Hd, see |8|. 
The Laplacian operator is discretised using the classical 7 point scheme, the discretised 
operator is designated in the following by A^. The total discretised magnetic field is 
then defined by 

H h (m) = AA h m + H d>h (m) + H a . h (m). 

Then, for a given external field I in Vl^R 3 ), we set nih,e, element of Wh(S 2 ), 
the equilibrium state of the discretised version of Q. This state is obtained using 
an explicit time discretisation combined with an optimisation of time which ensures 
its stability (see GOQOlEl). This equilibrium state, as seen previously, is such that: 

Vi G N h , 3/3, < Oand 

H h (m hA )li = f3imh,ili, 

we set Hh{mh,t) = Bg(mh,e) where Bi is a diagonal operator. Knowing an equi- 
librium state for the discretised sytem, we can define the linearised discrete system: 



(iuj - Dx th ,i{H h - Bt))nh = D\.h,iC,h (3) 

where D\ t h,e is the operator D\ t £ built for the mh,e equilibrium state. 

Then, for each element u of W^R 3 ), we associate a unique element U of R 3 #( Ar, >) 
defined by 



6 N h , Ui E M 3 and {/^ = — / u (x) dx. 

h 6 Jut 



Using this bijection between Wh(R 3 ) and R 3 ^^' 1 ), we can write a matricial version 
of the linearised discrete version of fl}: find U k in R 3 #( N ^ such that, for a given 
built on ^ we have 

M w U k =DY k , (4) 
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where, for every U in W?#( Nh \ for every i in Nh 



{MJJ)i 

and 




For use in the remainder of this paper for every U in R 3 #( JV ") we set: 

DHU = -MJJ + iwU 
then H is the matrix associated to the discrete operator Hh — Bg. 

2.3 Some properties of the discrete system © 

We set Mi, the element of R 3 *^) associated to mh,i- Let us define [m^] 1 " by 

W 1 = {w e c 3# ^)|Vz e iV fc , M^.Wi = o}, 

and we designate by the projection from £?#{ N h) m ( Then we can demon- 
strate: 

Theorem 1 For every Y in R 3 #( JVfc ) and for every lu strictly positive, the system (0 is 
regular and its solution is in an element of [rtih,e\ ■ 

Proof: If U is the solution of (|4}, then we have 

iuU = D(Y k +HU), 

knowing that D sends elements of C 3 #( N h) in [mf l .e]' L , we conclude that U is also an 
element of [mh.e] ± - 

Then, considering V in [m/^] , due to the structure of H, we have HV as an 
element of [m^/] . Each diagonal block (3x3) has 0,a + i and a — i as eigenval- 
ues. Knowing that the eigenvalues of H (symetric matrix) are real, we deduce that 
the eigenvalues of DH are complex numbers of non vanishing real parts unless the 
eigenvalueis null. Then, the eigenvalues of M u can not vanish. 

□ 

The conditioning number of the matrix M u can be estimated 
Theorem 2 For every u> real strictly positive, we have 

cond(M „) < ,/" 2 + " + ° 2 » 1 + *"'» + 1 + K » 
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Proof: This theorem is proved using the Courant-Fisher theorem for hermitian matrices 
which provides formulae for the highest and lowest eigenvalues. The proof is then 
classical and uses the fact that D*D is the projection matrix on [rat e] multiplied by 
(1 + a 2 ). 

□ 

We notice that the conditioning number cond(A/ w ) bahaves as expected when to tends 
to infinity: 

lim cond(M w ) = 1. 

ui — >oo 

Here, the fact that lo grows to infinity means that it dominates -p- . Now, if we consider 
that ijj is fixed, the behaviour of cond(M^) shows that the system is ill-conditioned 

lim cond(M w ) = oo. 

Thus, the pre-conditioning of the system is essential. In fact, the most interesting part 
of the spectrum of susceptibility for numerous applications is the low frequency part. 



3 The precontioning strategy 

3.1 Choice of the inversion method 

In order to solve system 0, we chose an iterative method; this choice is conditioned 
by the fact that the matrices considered are non-symmetric full matrices and the order 
of the systems to solve is great (up to 10 6 ). Three main iterative methods are used 
commonly to solve non symmetric systems: 

• the normal conjugate gradient (CNG), 

• the generalised minimal residual method (GMRES), 

• the conjugate gradient squared (CGS). 

As shown in the article of Nachtigal, Reddy and Trefethen II II . none of this three 
methods could be considered as a cure-all for all non-symmetric systems. As the con- 
vergence quality of CGS and GMRES is influenced by eignevalue clustering of the 
system matrix, CNG method convergence depends on singular value clustering. As the 
preconditionning strategy presented in this article is based upon the amelioration of the 
singular value clustering, we chose, of course, the CNG method. Furthermore, tests 
not presented in this article show that the CNG method seems to be more adaptated for 
this type of system, even if not preconditioned. 

3.2 An example of singular value repartition and of CNG conver- 
gence rate 

In the remainder of this paper, we have chosen to illustrate the results presented using 
a plain example. This example has been chosen for the low order, 192, of its system 
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which facilitates the visualisation (done with Matlab). The mesh chosen is a 4 x 4 x 4 
regular cubic mesh of a cubic domain. We set it in the dimensionless system A = 
0,88 1CT 10 , K = 0,57 10~ 2 , a = 0,5 and the cube length is equal to 10~ 6 . For 
this bench, we would want to choose to between io m i n — 0, 452 10 3 Hz and uj max = 
0,452 10 5 Hz. In Fig. |2]the error evolution for the CNG is shown. Here we have 

1 120, 1 1 1 1 1 1 1 1 — I 

bo 
O 

1 100- 




singular values 



Figure 1 : Singular value decompostion for the non-preconditioned system. 

chosen a final error criteria of 10~ 5 . With no preconditioning, the system converges in 
56 iterations for Lo m . m and 48 iterations for oj ma x, the precontioning number is almost 
equal to 7500 (slight variations between w m j„ and io max ). 




Iteration 

Fi gure 2: Residue of the CNG for 
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3.3 The preconditioning strategy 

We have three main goals to build the preconditioner: 

• to use the known properties of the system, 

• decrease the conditioning number sensitivity to the mesh size, 

• build a cheap preconditioner (memory size and computational cost). 

The first point is taken into account by using the result presented in Theorem ^ the 
right side of the preconditioner will be a projection on [mh,e) ± - This first step of 
projection eliminates the cluster of singular values near and ensures a convergence in 
48 iterations for uo min and of 27 iterations for Lo ma x- 



3.4 The "exact" preconditioner 

As a first stage, we would want to build a symmetric left precontioner. The non sym- 
metry of M u comes from the operator Di,h,e- In fact, we have 

D\,h,i w — —mi Aw + aPfW, 

the first part of the operator is a rotation, and the second part a projection. Then, it 
is possible to prove that D\,h,i does not have a main influence on the singular value 
decomposition. This means that we may choose a left preconditioner M g<u built on the 
operator 

iuj - a(H h - Be), 
That is to say, if we set H the matrix built on the operator Hh — Be 

Mg jUl = iuild — aH. 

In the sequel, even if we do not write the projection to lighten the notations, we consider 
that the system is right preconditioned by . Then, we prove the following theorem 

Theorem 3 For each uj and h strictly positive, we have 



condiM-lM^) < y/l+g(h,Lj) [2 + g{h,Luy 

where 

g(h,uj) - 



Proof: In the space [m/^] , we have: D = R + aid, where Id is the eye matrix on 
space [m/^] and R is the matrix associated to the operator — m^A. Moreover, for 
every U in [m,h.e] ± , we have 

(HU, m M ) = (U, Hm h . e ) = 0, 
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this implies, by breaking off of the elements of [mh t e}' L , that HU is an element of 
[mh,(\' L - We remark also that by working in [m,h,e]~ L , we have 

R 2 = -Id. 

Then, we have 

M-*M U = (iujH- 1 - Id)H~ 1 NH, 

where 

N h = -R H + (1 - a)H = NH. 
Then, for every V in [mh,e]' L , we have the following relation 

M-^MJl.M-^MJl = \\V\\ 2 + 2TZ[V.M-^N h V} + \\M-^N h V\\ 2 . 

We designate as 1Z[z] the real part of a complex number z. 
Furthermore, we have the following estimations: 

IKH^RHfW = UH^R 2 H) 2 \\ = 1, 
this implies that \\H~ 1 R H\\ = 1. So, using the fact than 

H~ 1 NH = -H~ 1 RH + (1 - a)Id, 

we have 

WH^N H\\ < + 

and 

L0 2 

\\M~IW = max(- T + l)- 1 <g(h,w). 
where Aj is the eigenvalues of the matrix H in [me^f l ] ± and is the set of indeces of 

Then, using the lowest eigenvalue controlled by projection part of the precondioner 
we conclude the proof of the Theorem. 

Finally, we have the good behaviour of the preconditionned system when the mesh 
length h tends to 0: 

lim (M~iM u ) < 1 + v/2 + a 2 , 

This version of the preconditioner gives excellent control of the conditionning number 
but needs the inversion of a full matrix. This leads us to the second stage in which we 
will replace the complete operator Hh by its laplacian part. 
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3.5 Preconditioning by the Laplacian component: 
proach 



the direct ap- 



The Laplacian part of Hh is the most punitive part of the matrix M w in terms of pre- 
conditioning. The idea in this section is to develop an approximate conditioner M g ^,A 
built on the operator 

ilu - AA h - B. 

The matrix M g<Ut j\ is a band matrix which could be more easily handled than M BtU , the 
earlier version of the preconditioned This approximation of the preconditioner M g>u 
will be all the more accurate as the norms of the operators Hd and H a are dominated by 
A. As seen in Fig. [3] the clustering of the singular value decomposition obtained for 
the system preconditioned by M gLU & is good. The convergence of the CNG algorithm 



5 1 20 




0.0625 



0.5 1 2 
singular values 



Figure 3: Singular value decomposition of the system preconditioned by M SiW ,a. 
using this method is very good: 7 iterations for u) m i n and 9 iteration for uj max (see Fig. 

n. 

3.6 The approximated preconditioner 

Nevertheless, the use of the pre-conditioner M g , u ,/\ stays expansive. The solution is 
to build an easily invertible approximation of M 9ia)i A- Here we will use here the work 
of 1 12 1. The idea is to project the matrix Mg.^./s into the circulant matrix space in the 
sense of the Froebenuis norm. 

Given a circulant matrix C n by n on C generated by c, vector of C n , we have 

V(z, j) G {1, n} 2 andp S {1, n}, Ci.j = c p ifj—i=p — l or j — i = n — pi. 

Then, as shown in 1 12 1, for every matrix M n by n on C, the projection C of M on the 
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Figure 4: Residue of the system preconditionned by M ffiW ,A- 



space of the circulant matrices of order n is generated by the vector c given by 

j / n—p+l n 

Vp € {1, ...,n}, Cp = - ^2 M l,l-p+l+ ^2 M l,l+n-p+l 
\ 1=1 l=n-p+2 

The three dimensional projection is more complex but the main idea is contained 
in the one-dimensional projection. 

When the circulant approximation matrix is built, the inversion is performed in 
the Fourier space (the matrix produced is block-diagonal 3 x 3 in Fourier space), then 
the precondioning is of complexity 0{N \og(N)) for each iteration of the inversion 
method. The other main advantage of the method is that the storage is reduced to 
O(N). 

In this section, we have to keep in mind that the structure is a three dimensional 
one: the considered matrices are 3 level block matrices. This implies that the projection 
must be performed on the 3 levels block circulant matrices. 

In the small example presented to illustrate the paper, the system precontioned by 
the approximated preconditioner converges in 30 iterations for the smaller frequency 
and 26 iterations for the highest (see Fig. |6j. The convergence curve is very good in 
the sense that the slope is quasi-constant. This point is quite important: susceptibility 
computations do not need high numerical accuracy. Effectively, the results obtained 
will be compared to experimental results for which the error is quite important. This 
comes form the fact that the samples used for experiments are far to be perfect and that 
the measurement tools do not have very high precision for this type of experiment. 
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Figure 5: Singular value decomposition of the system preconditioned by the circulant 
approximation of M 3yU) ,/\- 

4 Numerical simulations 

We present here the number of iterations for the simulation of a ferromagnetic dot. 
This dot is meshed by a regular grid, size 64x64x32. In this monolith a cylinder with 
the axis z and a circular basis (32 cells for the z direction and 64x64 for the others)is 
included. The total number of degrees of freedom is 393216. The results shown here 
have been computed on the parallel machines of ONERA and Dassault Aviation. 



4.1 Parallel implementation 

There are two possible levels of parallelisation for this problem: local parallelisation 
for computations of each iteration and global parallelisation of the frequency compu- 
tations. 

The global parallelisation is a repartition of each frequency computation through 
the processors. A main process distributes the computation to each processor such that 
each processor is always occupied. This part is implemented using MPI. 

The local implementation, not used for the results presented here, is the parallelisa- 
tion of the total magnetic field over the domain. In this computation, one part is more 
expensive than the others: the demagnetisation field. In fact, the computation of de- 
magnetisation is accelerated by using its Toeplitz structure (see 1 8 1). This computation 
strategy uses 3 dimensionnal FFT intensively. To enhance the performance, we have 
to parallelise the FFT computation. To do so, we have chosen to use OPEN-MR This 
choice avoids the transposition of the data via the cluster that must be performed while 
using a distributed memory system. The results are very satisfying: for a cubic struc- 
ture and sufficient number of cells (for instance a 32 x 32 x 32 mesh), the computation 
time of FFT is divided by 1.9 on a node of two processors. 
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Iteration 

Figure 6: Residue of the system preconditioned by the circulant approximation 

4.2 Description of the benchmark statistics and results 

The aim is to compute the susceptibility of a cylinder of permalloy (see for example 
1131 for this type of results). The parameters of the material are the following: 



Parameter 


Value 


A 

a 


0.17875 10" 11 
0.05 



In the following tabled we gi ye the number of iterations for directions x and y. 
The direction z in this computation is omitted because there is no resonnance in this 
direction. 

The computation has been carried out on a node composed of 8 Power4 IBM 
(1.1 GHz) with 16 GO of Ram. An iteration takes almost 24 seconds, the complete 
computation took 36 hours. 



w/1,356 (Hz) 


iterations for x 


error 


iterations i for y 


error 


3.00 10 b 


57 


4.92 10- 2 


127 


4.98 10~ 2 


2.73 10 5 


100 


4.85 10~ 2 


137 


4.91 10~ 2 


2.49 10 5 


198 


4.90 10~ 2 


267 


4.99 10~ 2 


2.26 10 5 


146 


4.99 10~ 2 


272 


4.95 10~ 2 


2.06 10 5 


188 


4.99 10~ 2 


329 


4.90 10~ 2 


1.88 10 5 


290 


4.96 10~ 2 


356 


4.91 10~ 2 


1.71 10 5 


316 


4.99 10~ 2 


317 


4.85 10~ 2 


1.55 10 5 


326 


4.98 10~ 2 


386 


4.94 10" 2 


1.41 10 5 


390 


4.99 10~ 2 


355 


4.94 10" 2 


1.29 10 5 


298 


5.00 10~ 2 


376 


4.97 10~ 2 


1.17 10 5 


329 


4.88 10~ 2 


354 


4.88 10" 2 
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iteration*; "for t 


error 


iteration*; i "for it 


error 


1 .U 1 _LU 


490 


4 97 1 0~ 2 

t-. y i _LU 


986 


4 91 1 0~ 2 

T^.7J -LU 


q 7 1 1 n 4 

y. i i j_u 


61 5 

ULJ 


4 91 1 O -2 

"-r.y i _LU 


400 


5 00 1 -2 

J.UU -LU 


s »4 1 n 4 

O.O^ -LU 


664 


4 99 1 -2 

"-r.yy -LU 


504 

Jut 


4 95 1 -2 

"-r.y^j xu 


s os 1 n 4 

O.UJ -LU 


M8 

U JO 


4 Q4 1 n -2 

l • /T- _LU 


41 6 
tiu 


4 811 n -2 

O 1 -LU 


7 w 1 n 4 

/ JJ -LU 


416 


4 Q4 1 n -2 

t."t ±u 


171 

j / 1 


4 Q9 1 n -2 

'-r.yz. ±u 


fi f>7 1 n 4 

u.u/ ±u 


118 

J 1 


4 Q6 1 n -2 

t.7U iu 


^1 Q 


4 70 1 n -2 


fi 07 1 n 4 


?Q1 


4 Q8 1 n -2 

T-.70 IU 


9Q4 


4 Q9 1 n -2 

T-.7Z ±U 


S 1 4 

J.JJ ±u 


1S1 


4 8S 1 n -2 

t.OJ ±u 


iUU 


4 qi i n -2 

'-r.y I ±u 


S 0^ 1 4 

-J.\JJ ±u 


177 

Oil 


4 Q1 1 -2 

t.yl -LU 


9S8 


4 78 1 n -2 

*+. / O -LU 


4 S8 1 4 

T.JO ±u 


411 


4 QQ 1 -2 
'-r.yy -LU 


9S9 


4 Qfi i n -2 

'-r.yxj ±u 


4. 1 7 1 4 
T. 1 / ±u 


tou 


4 Q7 1 -2 
'-r.y 1 _lu 


948 


4 8S 1 n -2 

t.OJ ±u 


^ 7Q 1 4 
j . i y ±u 


S41 


4 Q1 1 n -2 

'-r.y 1 _lu 


947 


4 7S 1 n -2 

'-r. I J _LU 


^ 4.S 1 4 


SQ9 


4 QS 1 -2 

t-.;7J -LU 


948 


s oo i n -2 

J.UU -LU 


^ 14. 1 n 4 


S47 


4 QQ 1 n -2 

'-r.yy _lu 


948 
z/+o 


4 qi i n -2 

t.7l ±u 


9 86 1 f) 4 


S4Q 


4 86 1 n -2 
t.ou ±u 


948 


s oo 1 n -2 

J.UU ±u 


2.61 10 4 


S71 


4 91 1 D -2 


247 


4 97 1 -2 

'-r.y i _LU 


1 37 1 4 

Z..J / _LU 


61 8 

UlO 


4 99 1 O -2 

'-r.yy iu 


?43 


4 99 1 n -2 

'-r.yy iu 


2.16 10 4 


653 


4 90 1 -2 


244 


4 99 1 -2 

'-r.yy iu 


1 97 1 4 

j. ,y i -LU 


686 

UOU 


4.62 10~ 2 


247 


4 99 1 n -2 


i 7Q i n 4 

i . / y _lu 


791 

/ z, j 


4 8^ i n -2 

T.OJ L.U 


?S1 


4 Q9 1 n -2 

'-r.yz. ±u 


1 m in 4 

l.UJ -LU 


77Q 

i i y 


4 QO 1 n -2 

T.J/U l.u 


?S4 


4 Q7 1 n -2 

*+. y I ±u 


i 48 in 4 


81Q 

ojy 


4 Q4 1 n -2 

L +.y-T l.u 


?S7 


4 QQ 1 n -2 

'-r.yy ±u 


1 is in 4 

1 . J J _LU 


8SS 

OJJ 


4 QO 1 n -2 

T.J/U L.U 




4 86 1 n -2 
t.ou ±u 


1 ?i in 4 

1 .Z. J J.U 


842 


4 99 1 -2 

"-r.yy iu 


267 


4 9S 1 -2 

t.7J -LU 


1 191 n 4 


819 

O JZ. 


4 QO 1 n -2 

^■."W -LU 


97^ 


4 QS 1 n -2 

t.7J ±U 


1 n? in 4 

1 .uz. _lu 


819 


4 Q4 1 n -2 

t."t ±u 


97Q 


4 Q6 1 n -2 

t.7U ±U 


q ii 1 n 3 

y .z. i _LU 


81S 


4 Q9 1 n -2 

'-r.yz. _lu 


980 


4 Q8 1 n -2 

t.7u ±u 


8 44 1 n 3 

O .ii _lu 


6Q1 

U~ 1 


4 Q9 1 n -2 

'-r.yz. _lu 


98Q 


4 8Q 1 n -2 

'-r.oy ±u 


7.68 10 3 


843 


4.96 10~ 2 


292 


4.94 10~ 2 


6.99 10 3 


856 


4.98 10~ 2 


296 


4.97 10" 2 


6.36 10 3 


868 


4.90 10~ 2 


302 


5.00 10~ 2 


5.79 10 3 


875 


4.93 10~ 2 


305 


4.94 10" 2 


5.27 10 3 


888 


4.83 10~ 2 


311 


4.95 10" 2 


4.80 10 3 


893 


4.98 10~ 2 


316 


4.99 10~ 2 


4.37 10 3 


907 


4.94 10~ 2 


319 


4.94 10" 2 


3.98 10 3 


913 


4.92 10~ 2 


327 


4.94 10~ 2 


3.62 10 3 


923 


4.97 10~ 2 


328 


4.99 10" 2 


3.30 10 3 


944 


4.99 10~ 2 


336 


4.91 10" 2 


3.00 10 3 


951 


4.85 10~ 2 


339 


4.94 10~ 2 



Tabled Iteration table. 
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5 Conclusion 



The goal of the study was to allow the computing of micro-wave susceptibility of ferro- 
magnetic particles with thin details. This last point called for very large meshes (about 
300000 degrees of freedom) for which the classical inversion methods with no precon- 
ditionning did not work at all, or required such a large amount of iterations that the 
computation times for an acceptable range of frequencies was far from useful. The 
strategy presented in this article is an industrial computations approach, and obtains 
interesting results for a large spectrum of benchmark. Computations of realistic ex- 
periments have been performed (see I14II13|[T31 ) for physical systems where it was 
possible to compare results with physical experiments. Some problems remain, in par- 
ticular, the strategy developed aims at the laplacian part of the total magnetic field 
whereas some systems are revealed to be principally influenced by the demagnetising 
field. The next step is to extend the strategy of the paper in order to include the de- 
magnetisation part of the magnetic field in the approximated preconditioner. The main 
difficulty of the extension is algorithmic: to build a good circulant approximation of 
block Toeplitz matrices. An another interresting point to study would be the implemen- 
tation of an efficient parallelised FFT algorithm for distributed memory systems. The 
main problem of such an implementation would be the optimisation of the transposition 
phase of the data through the memory nodes of the distributed system. 
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2.5 



Relaxation pour un preconditionnement circulant sans projection 




Numero d'iteration 



Decomposition en valeur singulieres de la matrice non projetee preconditionnee 



frequence la plus haute 
frequence la plus basse 
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Decomposition en valeur singulieres de la matrice projetee preconditionnee 



frequence la plus haute 
frequence la plus basse 
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Decomposition en valeur singulieres de la matrice projetee preconditionnee 
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